Enzymatic Hydrolysis of Human Milk Oligosaccharides. The Molecular Mechanism of Bifidobacterium Bifidum Lacto-N-biosidase

Bifidobacterium bifidum lacto-N-biosidase (LnbB) is a critical enzyme for the degradation of human milk oligosaccharides in the gut microbiota of breast-fed infants. Guided by recent crystal structures, we unveil its molecular mechanism of catalysis using QM/MM metadynamics. We show that the oligosaccharide substrate follows 1S3/1,4B → [4E]‡ → 4C1/4H5 and 4C1/4H5 → [4E/4H5]‡ → 1,4B conformational itineraries for the two successive reaction steps, with reaction free energy barriers in agreement with experiments. The simulations also identify a critical histidine (His263) that switches between two orientations to modulate the pKa of the acid/base residue, facilitating catalysis. The reaction intermediate of LnbB is best depicted as an oxazolinium ion, with a minor population of neutral oxazoline. The present study sheds light on the processing of oligosaccharides of the early life microbiota and will be useful for the engineering of LnbB and similar glycosidases for biocatalysis.


System preparation
The initial coordinates for the simulations were taken from the crystal structure of Bifidobacterium bifidum's LnbB in complex with LNB-thiazoline (PDB accession code 4JAW, resolution 1.80 Å) 1 (Figure 1). The model was built using chain A. Residues 30 to 40, which they were added in the cloning process for a better packing of the construct, were removed, as well as the LNB-thiazoline ligand. The protonation states of aspartate, glutamate and histidine residues was assigned considering the environment of each residue and the optimal pH for the enzyme activity (pH 4.5) which was also tested with PropKa3 and H ++ 2-4 . The assisting residue (Asp320) was considered as deprotonated, and the acid/base residue (Glu321) was taken in its protonated form. The LNT substrate was docked in the active site using the program AutoDock VINA 5 ( Figure S2). Different poses were generated and the best ones in terms of contacts with the residues of the -2 and -1 subsites (reproducing the ones observed in the X-ray structure of the enzyme with LNB-thiazoline, PDB 4JAW) were kept and submitted to classical MD simulations.

Classical molecular dynamics
The AMBER18 6 package was used in all classical molecular dynamics simulations. The system was solvated in a cubic cell of 33868 water molecules and 10 Clanions were added to neutralize the system. The force fields FF99SB 7 , GLYCAM_06j 8 , and TIP3P 9 were used for the protein residues, the oligosaccharide substrate and the water solvent molecules, respectively.
Classical MD simulations started with energy minimization. First, only water molecules and Cl − anions were minimized (from residues 629 to 34507), followed by minimization of the whole system, using steepest descent and conjugate gradients methods. Next, the system was gradually heated to 300 K by increasing the temperature 100 K over 50 ps runs (only solvent and ions during the first 100 K, then the whole system until 100K, 200K and 300K), using a timestep of 1 fs. Once the system was at 300 K, the density was equilibrated in the NPT ensemble (2 ps). Finally, the simulation was extended in the NVT ensemble for 50 ns ( Figure   S3 S3) using a timestep of 2 fs. Restraints at interaction distances in the active site were used to stabilize the ligand at the binding site, which were gently released during the production run.

QM/MM molecular dynamics
QM/MM MD simulations were performed using the CPMD code, version 3.15.1 (http://www.cpmd.org). 10 A snapshot of the classical MD of the enzyme-substrate complex was used as starting point for QM/MM. The QM region was enclosed in a 15.09 × 19.84 × 18.52 Å cell with the side chains of the catalytic residues Asp320 and Glu321, the entire GlcNAc unit at the -1 subsite and part of the sugars at the -2 and +1 subsites, in total 67 atoms including 6 monovalent link atoms that were used to saturate the QM region 11 ( Figure   S5). The remaining 111345 atoms were treated with molecular mechanics (MM). A larger QM region including the positively charged His263 residue was also tested to make sure that there is no proton transfer between this residue and Asp320 in the Michaelis complex.
QM/MM simulations were performed using the method developed by Laio and coworkers, 12 which combines Car-Parrinello MD with DFT. 13 The QM/MM interface was treated using the multilayer electrostatic coupling scheme developed by Laio and coworkers, 12 with NN, MIX and ESP radii of 9.52, 11.11 and 15.87 Å, respectively. 11 We used the Perdew-Burke-Ernzerhof (PBE) generalized gradient-corrected approximation for DFT 14 , as it has been assessed a suitable option for describing carbohydrate conformations 15 and other glycoside hydrolases. [16][17][18][19] A 70 Ry kinetic energy cutoff we used for expanding Kohn-Sham orbitals in the plane-wave basis set, along with norm-conserving ab initio Trouiller-Martins pseudopotentials. 20 The fictitious electron mass was set at 600 au and the timestep at 0.12 fs.
Geometry optimization was carried out with QM/MM with annealing of the ionic velocities until the nuclear gradient is lower than 5·10 -4 au. The Nosé Hoover thermostat 21 was used to heat the system to 300 K during a 5 ps MD run. The final snap-shot of the QM/MM MD simulation was taken to start metadynamics simulations.

Classical molecular dynamics
The enzyme-intermediate complex was solvated in a cubic cell of 33883 water molecules and 9 Clanions. AMBER 18 was used in all simulations, along with forcefields FF14SB 28 (protein) GLYCAM_06j (for the Gal part of LNB-oxazoline) and TIP3P (solvent).
Parameters were assigned according to the GLYCAM_06j and GAFF 30 force fields, using the Antechamber suite 31 of AMBER18. RESP charges for the LNB-oxazolinium were obtained using Gaussian09 29 at the HF/6-31G* level of theory Modeling of the reaction intermediate with classical MD was performed using the protocol described in section 1.2. A detailed analysis of the dynamics of the catalytic water molecules moving near the reactive center during the MD simulations was performed using an in-house Python3 program. To identify water molecules well-positioned for catalysis, a cutoff of 4 Å was chosen for the water nucleophilic attack distance at the C1 and 2.2 Å for the hydrogen bond with the general base Glu321. One of the structures fulfilling these criteria was chosen as starting point for the QM/MM calculations.

QM/MM molecular dynamics
The QM region (56 atoms) was selected as including the side chains of Asp320, Glu321, the catalytic water, the GlcN-ox + sugar ring at the -1 subsite and part of the Gal sugar at the -2 subsite ( Figure S5). The QM/MM boundary was treated with monovalent link atoms as previously described. The QM region was enclosed in a 16.50 x 19.97 x 15.39 Å cubic cell.
The MM region contained 111354 atoms, which were treated classically. The protocol used for QM/MM MD simulations was the same as described above (section 1.3). The radii of the zones of electrostatic coupling was taken as 15.34 Å (ESP), 11.64 Å (MIX) and 8.46 Å (NN).
The fictitious electron mass was set at 600 au and the timestep at 0.12 fs. A geometry optimization was carried out by annealing on the ionic velocities until the nuclear gradient is lower than 5·10 -4 au. The system was heated to 300 K by coupling it to a Nosé-Hoover thermostat for 4 ps.

QM/MM metadynamics
QM/MM metadynamics were used to model the second reaction step leading to product formation. CVs including broken and formed bonds were used. One CV was taken as

Cartesian coordinates of the main states along the reaction pathway
The Cartesian coordinates (PDB format) of the atoms included in the QM region in states MC', MC, TS1, INT and P (representative structures shown in Figure 3C) are listed below.